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Abstract 



We show that the flux- field expansion derived by Boghosian & Covcney ( 200Cl| ) for the 



Rothman-Keller immiscible fluid model can be derived in a simpler and more general way in 
terms of the completely symmetric tensor kernels introduced by those authors. Using this gen- 



eralised flux-field expansion we show that the more complex amphiphilic model of Boghosian 



Covcney & Emcrton (1996) can also be derived from an underlying model of particle interac- 
tions. The consequences of this derivation are discussed in the context of previous equilibrium 
Ising-like lattice models and other non-equilibrium mesoscale models. 
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1 Introduction 

Modelling and simulation of amphiphilic fluids remains a challenging area of physical and com- 
putational science. In equilibrium, ternary amphiphilic fluids self-assemble into a multitude of 
mesophases. The morphology of each mesophase depends on the molecular character of the am- 
phiphile as well as on temperature, concentration, and solvent salinity. In turn, the mesophase 
morphology influences the nonequilibrium properties of the fluid. Wormlike micellar and sponge 
phases may exhibit visco-elastic rheological properties, shear induced phase transitions or demix- 
ing under flow. For a comprehensive review of amphiphilic behaviour and (mostly equilibrium) 



modelling techniques see [Gompper fc Schicl^ (|1994D . 



These examples illustrate the central challenge posed by the study of amphiphilic fluids: how 
does one construct models which are sufficiently computationally tractable to reach the time and 
length scales of interest (/im - cm) while including sufficient detail from the molecular scale (nm) 
to describe the behaviour of specific compounds? 

The present article goes some way to addressing this issue. In particular, it shows how a previ- 
ously phenomenological lattice-gas model of amphiphile dynamics can be related to an underlying 
model of particle interactions. In section |2| we describe previous hydrodynamic lattice-gas automata 
models used to simulate binary immiscible and ternary amphiphilic behaviour. In section ^ we de- 
scribe the particular collision rules used in the Boghosian, Coveney and Emerton (1996) model. 
In section H we specify the form of the underlying model of particle interactions and obtain an 
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expression for the change in system energy due to both colhsion and propagation. In section |5| we 
obtain the generahsed form of the Boghosian &: Coveney| ( 2000| ) flux- field expansion and use it to 
show that the Boghosian, Coveney and Emerton (BCE) model is a truncation of this expansion at 
first order. 



2 Hydrodynamic Lattice-Gas Automata 



Lattice-gas automata have been used extensively for modelling hydrodynamics since Frisch, Has- 
slacher, and Pome au ([Frisch, d'Humieres, Hasslacher, Lallemand, Pomeau &: Rivet 1987] ), and 
Wolfram ( Wolfram 1986| ) showed that it is possible to simulate the incompressible Navier-Stokes 
equations using discrete Boolean elements on a lattice. The dynamics of all hydrodynamic lattice 
gases take place in two substeps: all the particles simultaneously propagate along their velocity 
vectors to arrive at a new lattice site (retaining their velocity vectors as they do so), and then col- 
lide with the other arriving particles. The collisions are required to conserve mass and momentum. 
Even under the constraints imposed by conservation laws, a multiplicity of collision outcomes are 
usually possible. 

The simplest adaptation of the single phase lattice-gas which enables one to simulate multiphase 
phenomena is the introduction of an internal degree of freedom (termed 'colour') for the particles. 
For example, if the particle colour takes values (-1-1,— 1), one has two phases, which are typically 
referred to as red and blue. In order to extend the range of multiphase behaviour captured by the 
lattice-gas, it is necessary to include colour dependent collisions. This was first done by Rothman 
and Keller in order to allow simulation of immiscible fluids (Rothman & Keller 1988| ). Rothman and 



Keller introduced a colour field, which gives information about the composition of the neighbour- 
hood of a site, and a colour flux. Those authors introduced a simple collision rule, choosing collision 
outcomes based on the colour flux and colour field which have the effect of inducing cohesion in 
the red and blue phases. 

The Rothman and Keller model has been used to investigate surface tension and interfa- 



cial fluctuations ( Adler, d'Humieres &: Rothman 1994 ), spinodal decomposition in both two and 
three dimensions (Rothman & Keller 1988, Appert, Olson, Rothman &; Zaleski, 19951 ), and the 
effect of shear on phase separation in three dimensions ( Olson &: Rothman. 1995| ). [Emerton " 



Coveney &: Boghosian ( 1997g| ) generalised the Rothman-Keller model by including point surfac- 
tant particles, possessing a colour dipole. Those authors introduced three additional terms into 
the Rothman-Keller model. The form of these terms was based on the dipole-charge, charge- 
dipole and dipole-dipole interactions. The terms were then combined in a local Hamiltonian 
with arbitrary parameters, which were subsequently fixed by an extensive parameter search for 
canonical amphiphilic behaviour (micellisation) . The derivation of their local amphiphilic Hamil- 
tonian is described in ( Boghosian et aT] [1996[ ) and its form is given below. This model has 
been extensively studied in both two and three dimensions, and has been very successful in cap- 



turing a wide range of both equilibrium and non-equilibrium surfactant behaviour ( Emerton et 
n [1997c, Coveney, Emerton & Boghosian 1996, Emerton, Coveney & Boghosian 19976[, [Emerton 



oveney Sz Boghosian[ [19974 [Emerton, Weig, Coveney fc Boghosian[ [19974 [Boghosian, Coveney 



Love 2000, Boghosian, Coveney, Love fc Maillet[ [2001 , Love, Coveney &: Boghosian 20010, Love, 
aillet &: Coveney 2001J, Love &: Coveney 200 1[ ). However, the model lacks any clear connection 



to a molecular description of surfactants, and as such the connection between different parame- 
terisations of the model and different surfactant compounds remains unclear. Additionally, there 
exist many equilibrium lattice models (see pompper &: Schick[ ( [1994[ ) for an overview) which are 
superficially very similar to the Boghosian, Coveney and Emerton (BCE) model, and it is unclear 
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how the equihbrium phase diagram of the lattice gas model is related to the equilibrium phase 
diagram of those models. 

Recently Boghosian and Coveney derived the Rothman-Keller model from an underlying model 
of particle interactions ( [Boghosian fc Coveney) |2000 ). In the limit where the ratio of the mean 
free path to the interaction range (denoted by e) is small, those authors showed that an arbitrary 
potential may be expanded in terms of fluxes and fields. The Rothman Keller model is a truncation 
of such an expansion at first order in e. This work enables one to reinterpret previous phenomeno- 
logical lattice gas models in terms of an underlying potential. This has two consequences. Firstly, 
the interaction potential will have a specific range, which in turn associates a characteristic scale 
for one lattice spacing. Secondly, the parameterisation of the interaction potential enables one 
in principle to properly introduce molecular specificity. However, the derivation of the flux-field 
expansion from the underlying potential obtained for the Rothman and Keller model was of con- 
siderable complexity, making the direct extension of this work to more complex local Hamiltonians 
daunting. In fact, those authors noted that the simplicity of the final result obtained implied that 
a simpler derivation might exist. 

In this paper we show that there is indeed a simpler derivation of Boghosian and Coveney's 
result. We write our initial potential in terms of the symmetric tensor kernels introduced by 
Boghosian and Coveney. By directly considering the Taylor expansion of the mth rank symmetric 
tensor kernel, and utilising a recursion relation between the tensor kernels, we show that the flux- 
field expansion may be obtained much more directly. In addition, this derivation allows us to 
naturally generalise to the higher order tensor fluxes required for the amphiphilic lattice gas. The 
local amphiphilic Hamiltonian of the BCE model is obtained from this flux-field expansion to first 
order, in exactly the same way as the Rothman-Keller model was obtained for the binary fluid case. 



3 Lattice-gas amphiphilic fluid dynamics 

In this section we outline the original phenomenological BCE model. In single-phase lattice-gas 
dynamics a given incoming state s has a set of possible outgoing states {s'} which conserve mass 
and momentum. This set is referred to as the equivalence class of the state s. The hydrodynamics 
will be identical for any choice of s' from {s'}, but in a multiphase model, the colour dynamics will 
be different. In the Rothman-Keller model, for each member of {s'}, a colour flux is computed, 
and the collision process is modified such that the notional colour work done by the flux against 
the field is minimised. Chan and Liang introduced a sampling procedure over the equivalence class 
members such that s' is sampled from a probability density 'P(s'), which is modelled as the Gibbsian 
equilibrium corresponding to a Hamiltonian H{s'): 

V{s') = ^exp[-PH{s')], (1) 

where P is an inverse temperature, H{s') is the energy associated with collision outcome s', and Z 
is the equivalence-class partition function. Whereas the Rothman-Keller model had a single term in 
the Hamiltonian, the BCE model has three additional terms, capturing the interaction of surfactant 
particles at a site with neighbouring colour charges, the interaction of colour charges at a site with 
neighbouring surfactant particles, and the interaction of surfactants with neighbouring surfactants. 
It is these four terms in the Hamiltonian which we systematically obtain in the remainder of the 
paper. 
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4 Collisional Energetics 



We denote the postcollision charge-like and dipolar attribute with velocity Cj at site x by Q^(x) and 
p^(x). Upon subsequent propagation, the charges and dipoles (/^(x), Pi(x) will be at position x + Cj 
and the charge qj{x + y) and dipole Pj(x + y) will be at position x + y + Cj. This is illustrated 
in Fig. 1^. We take as our starting point the BCE potential energy for a set of point charges and 
point dipoles: 



+ 



\ [%(x)Pj(x + y) •y(/>i(|y| 
1 



9i(x + y)Pi(x) -y^idyl 



[pj(x)pj(x + y) :yy(/>2(|y|) + Pj(x)pj(x + y) : M\i\y\ 



Where we have defined the following functions related to the derivatives of (/>(y): 

<Am(y) ■ 



ydy 



(2) 



(3) 



In this case there are three contributions to the potential energy: the first term is the charge-charge 
interaction of the RK model, the second is the charge dipole interaction and the third term is the 
dipolar-dipolar interaction. We have explicitly separated the charge-dipole contribution into two 
terms, to make it clear that vector i at site x may be occupied by either a charge or a dipole. The 
factor of 1/2 premultiplying this term compensates for double counting. We now introduce the 
completely symmetric mth rank kernel defined by Boghosian and Coveney: 

^2m— n 



E 



(n — m) 



;per 




(4) 



m=[n/2] 

where "per" indicates a summation over all distinct permutations of indices. This enables us to 
write the BCE potential energy in a particularly compact form: 

+ fe(x + y)pi(x) - gi(x)pj(x + y)] •/Ci(y) 
- Pi(x)pj(x-hy) : /C2(y) 



(5) 



The change in the potential energy due to both collision and propagation is then given by: 



tot 



AV„. 



+ 



AVn + AVc 

^EENwM(^ + y)-'?^ 

A[p;.(x)p;.(x + y)] :/C2(y) 
i^^g:(x)g;(x + y)A/Co(y) 

[5i(x)Pj(x 



y)Ap^(x) ./Ci(y) 



+ y)-gUx + y)p:(x) •A/Ci(y) 



+ p;.(x)p;.(x + y) : A/C2(y), 



(6) 
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X + y + Cj 



X + y 




Figure 1: Change in Potential Energy: The change in the potential energy of interaction be- 
tween two particles moving in (possibly) different directions at (possibly) different sites is illustrated 
here. See Eq. (|6|). 

where we have separated the contribution to Ay due to the movement of the interacting particles, 
AVc, and that due to nonconservation of the dipole direction in collisions, AVn- 

5 The Flux-Field Decomposition 

Boghosian and Coveney proceeded by expanding A/Co(y) in a Taylor series in the ratio of the 
characteristic lattice spacing c to the characteristic interaction range y. We wish to expand A/Co(y), 
A/Ci(y) and A/C2(y) in a similar way. These first and second rank tensors are functions of the 
vector y, not only of its magnitude. The Taylor expansion of the mth rank kernel is given by: 

OO -. 

This may be written as a sum over nth rank inner products of nth rank outer products: 

OO , n \ n n 

^m(y + c) = ^ - (g) V /C,„(y) O (8) c, (8) 

On . \ / 

enabling one to use the property of the nth tensor derivative of the mth rank kernel (which is 
proved in Appendix^): 

((8)v)/C„,(y) =/C™+„(y). (9) 

Letting c = Cj — Cj in equation (^): 

OO 

lCm{y + Cj-Ci) = ^_[(c,--c,)-Vr/C™(y) (10) 

^ Ti. 
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and binomially expanding [{cj — Cj) • V]": 

[{cj - • V]" = X: {fj (-l)'(c, • Vr-Hc, • V)'}Cm{y), (11) 

writing (cj • V)"'"'^ and (cj • V)'^ in terms of nth rank inner products of nth rank outer products as 
in equation gives: 

OD -, n / \ / k \ ^ n—k /n~k \ 

/C^(y + c,-c,) = Y.^-^j:J(lY-ify^cj^ (12) 



Using (|l2|) to substitute for A/Cm,(y) in ^Vc give the charge-charge interaction as: 

A^ce = jEE3E^.^WOU M)'^EMy)0^n%(x + y) 

2 X n=l fc=0 V*^/ y 



and the charge-dipole interaction as: 

^^c, = ^EE^E^^'WO ! M)'E^"+i(y) O ^n%+i(x + y) 

2 X n=l^-k=Q \'^J y 

^ oo 1 n A;+l / \ n—k 

- i E E ^ E ^r^i(x) O U (-1)' ^ O ^n'-.(x + y) 

2 X „=i"-fc=o y 



and the dipole-dipole interaction as: 

oo 1 ?^ fc+l / \ n— fc+1 

^y^^ = iEE;^E^^+lOU M)'E^n+2(y) O ^n%+l> 

2 X „,=i"-fc=o VV y 



where we have defined: 



^i(x) = Ep^ 



k 

x)(8)Ci 



J,«(x) = 5:gKx)(g)c. 



(13) 



(14) 



(15) 



(16) 



Because these equations are summed over the whole lattice we have a freedom of choice in labelling 
the particles. This means that the above expressions are invariant under the interchange: 

i ^ j 

j ^ i 

X ^ X + y 

y ^ -y (17) 
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In the first and third terms this means that the kth term and the n — kth term are equal. For the 
charge dipole terms the kth term in the qipj series is equal to the n — kth term in the qjPi series. 
We can therefore drop half of the sum over k. To first order in epsilon the change in potential 
energy due to the motion of the particles is therefore: 

AyW = -^^''(x)Q5]/Ci(y)Jo^(x + y) 

y 

- :^/(x)OE^2(y)6jf(x + y) 

y 

+ :^|(x)6E^2(y)Jo'(x + y) 

y 

- ^l6E^3(y)6^f- (18) 

y 

Introducing the definitions of colour flux: 

J = J?, (19) 

colour field: 

E = -E^i(y)^o(x + y), (20) 
y 

dipolar field: 

P = -E^2(y)- Jflx + y), (21) 
y 

dipolar flux tensor: 

^ = ^2'' = EPi(^)^- (22) 

i 

colour field gradient tensor: 

^ = E^2(y)Jo'(x + y), (23) 
y 

and dipolar field gradient tensor: 

p = -E^3(y)- JiV + y)- (24) 

y 

from the BCE model. Substituting these definitions into the first order expansion (|l8|) yields: 

Ayj^) = 3 ■(E + P)+J -.{S + V) (25) 

which is identically the BCE local Hamiltonian without the terms due to collisional non-conservation 
of dipolar moment. 

It only remains to evaluate AV^. Firstly we retain only those parts of AVn which distinguish 
between outgoing states, writing: 

= EE'?i(^ + y)P^(^) • ^i(y) + I [p;-(x)p;(x + y)] : /C2(y) (26) 
x,y i,j 
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Using the definitions (18-23) and defining the total outgoing dipole vector: 

p'(x) = ^pKx) (27) 

i 

we obtain: 

A^n = ^Ep'W-(2E + P) (28) 

X 

This has exactly the form used in the BCE model. The factor of two above, which was not included 
in subsequent utilisations of the model, arises because they derived the various terms in a local 
Hamiltonian piecewise, without taking account of double counting considerations. These terms in 
the local Hamiltonian were then combined with arbitrary coefficients. 



6 Discussion and Conclusions 



Having obtained the BCE model Hamiltonian from an underlying potential, it is now possible 
to obtain parameterisations of this model for a given form of interaction potential. Equally, it 
is possible to attempt the inverse problem and use the much studied two- and three-dimensional 
parameterisations of the model to obtain a form for the interaction potential. Any computational 
implementation of the model will utilise particular stencils to evaluate the tensor fields and gradients 
over a finite range. These stencils, together with the form of the Hamiltonian defined here and a 
given model parameterisation will yield a set of differential equations and boundary conditions for 
(piy). Once obtained, this form of (/>(y) could be utilised in other mesoscale models. Most notably, 
the amphiphilic Malevanets-Kapral described in (Sakai, Chen & Ohashi 20001 ) is very similar to 
the BCE model. This raises the interesting possibility of studying the same underlying model of 
particle interactions in both lattice and off-lattice models. However, previous studies of the BCE 
used nearest neighbour stencils. The interaction range is therefore comparable to the mean- free 
path (epsilon is close to one) and so additional studies of the BCE model with longer-range stencils 
may be necessary to make this comparison. 

Additionally, many equilibrium lattice models of amphiphilic systems exist, whose phase di- 
agrams have been obtained by Monte-Carlo simulation. One in particular, that of Matsen and 
Sullivan utilises the same form for the potential of a set of interacting charges and dipoles as that 
used here (Matsen & Sullivan 199C, [Matsen, Schick &: Sullivan 1993| ). A given set of parameters 



for this equilibrium model can therefore be directly translated into a set of parameters for our 
non-equilibrium model. This enables the equilibrium phase diagram of a given set of parameters 
to be calculated by Monte-Carlo simulation, or equivalently the dynamic behaviour calculated for 
a set of parameters for which the equilibrium phase diagram is already known. 

By giving previously phenomenological models an underlying theoretical foundation, the con- 
ceptual problems inherent in the lattice-gas approach, and in mesoscale modelling in general, are 
thrown into sharp relief. Although we have referred to a 'local Hamiltonian', in fact the lattice 
gas algorithm does not pay any attention to energy conservation. One could consider following the 
same approach as that used for energy-conserving DPD, and attach an internal degree of freedom 



to the lattice gas particles (Espanol 1997). However, such an internal energy changes the nature of 
the particles from microscopic to thermodynamic entities. Their dynamics should then be required 
not only to satisfy energy conservation, but also guarantee positive entropy production. 



The lattice-gas models studied in (Emerton et al. 1997q, Coveney et al. 1996, Emerton et al. 


1997?^ 


Emerton et al. 


1997c, 


Emerton et al, 19974 Boghosian et al. 2000, Boghosian et al. 2001, Love 


et al. 


2001Q 


, Love et al. 


20016, 


Love & Coveneyl 2001|) have elements of both interpretations. On 
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the one hand lattice gas surfactant particles are regarded as being genuinely molecular in character, 
whereas in the collision step outgoing states are sampled over Boltzmann weights constructed from 
the local Hamiltonian. This sampling procedure effectively puts each lattice site in contact with a 
heat bath, thereby introducing a thermodynamic aspect to the simulation. This thermodynamics is 
flawed in three ways. Firstly, as was shown by Chan and Liang, there is a net heat flow between the 
heat bath and the simulation, despite the fact that they are at the same notional temperature ( phan 



fc Liang 1990 ). Secondly, the total energy of the heat bath plus system is not conserved, and thirdly 
the system does not possess a H-function. Put more simply, the thermodynamic aspects of these 
lattice-gas simulations violate the zeroth, first and second laws of thermodynamics. 

Thermodynamically consistent mesoscale models remain elusive. At present only the GENERIC 
formulation of DPD can claim complete thermodynamic consistency, at the expense of abandoning 



any molecular level information (Espanol, Serrano & Ottinger 19991) . Such a completely top-down 



approach can only be useful for materials for which a widely accepted macroscopic description 
exists. Such materials include liquid-gas mixtures or nematic liquids, but no commonly agreed 
upon continuum description exists for many materials, including amphiphiles. 

The thermodynamic difficulties of the lattice-gas algorithm described above arise partly from 
the mixture of microcanonical and canonical elements. In the context of equilibrium Monte-Carlo 
simulation an approach by Creutz shows how to perform equilibrium Monte-Carlo simulations in 
which microcanonical and canonical elements are incorporated in a systematic way (Creut2 1983). 



It is possible that similar approaches to the non-equilibrium simulations considered here could also 
yield useful modifications. 

Acknowledgements 

The author would like to thank Peter Coveney, Bruce Boghosian and Jonathan Chin for productive 
discussions and CINECA and the EU human mobility program for a MINOS grant during which 



part of this work was undertaken. Figure 1 is reproduced with permission from Boghosian & 
Coveneyl (^00C| ). 



9 



A Proof of recursion relation eqn (|^). 



We now prove the recursion relation @ by induction. We start by noting that @ is triviahy true 
for n = l,m = and therefore the result V/C„,(y) = lCn+i{y) completes our proof of (|9|). Writing 
out the left hand side of this expression and using V4>m{y) = y4>m+i{y) gives: 



V/C„(y) 



E 



m=\n/2\ 



4>m+i{y) 

(n — m)! 



yper [2m — n,n — m] 



+ E 



(n — m) 



■Vper [2m — n,n — m] 



m=\n/2\ 

where we have introduced the further shorthand: 

/2m— n 

per 



per [2m — n,n — m] . 



(29) 



(30) 



By substituting m' = jn — 1 in the first term here and considering series for odd and even n 
separately, we obtain: 



V/C„(y) 



+ 



E 



(n/2)! 

4>m{y) 



(n + 1 — m)! 



Vper [0, m] 

yper [2m — 2 — n, n + 1 — m] 



m=[n/2+l] 

+ (?i + 1 — m)Vper [2m — n,n — m] 



+ 



(pn + l(y)yper [n,0] . 



(31) 



for even n and: 



V/C„(y) 



E 



r{y) 



(n + 1 — m)] 



yper [2m — 2 — n, n + 1 — m] 



m=\n/2+l] 

+ (?i + 1 — m)'Vper [2m — n,n — m^ 



+ 



<Pn + l{y)yper [?i,0] . 



(32) 



for odd n. Using yper [n, 0] = per [n + 1,0] and Wper [0, m] = enables us to reunite these two 
series, giving: 

V/Cn(y) = V 7 — ^"^^^'^ yper [2m + 2 - n,n + 1 - m] 

m=r(n+l)/2] ^ 

+ - — ^'^^y^ — + 1 — m)Vper [2m — n, n — ml 

+ (/.,+i(y)per[n + l,0]. (33) 
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We can simplify this expression if we consider the types of additional permutations of indices in 
per [2m — n — l,n + 1 — m] as compared with those in per [2m — n,n — m]. Clearly, the (n + l)th 
rank tensor has one additional index, which can be attached either to a y or to a Kronecker delta. 
If the new index is attached to a y there is only one additional permutation, which is taken care of 
by the first term in (|3^). The second term will result in 2m — n terms on differentiation, the new 
index being combined in each successive term in a Kronecker delta with each of the 2m — n indices 
labelling the y's in the original expression. These terms capture all the unique ways of assigning 
the new index to a Kronecker delta. 

The prefactor (n + 1 — m) is necessary as the whole expression is divided by {n + 1 — m)!, 
which arises as this factor is exactly the number of ways of permuting n+1 — m Kronecker deltas, 
this corresponds to (n + 1 — m)\ equivalent permutations of the indices. These permutations are 
not generated by the outer derivative, but this factor captures their effect on the resulting tensor, 
enabling us to write: 

i^per [2m — 2 + n, n + 1 — m] + Vper [2m — n,n — m]^ = per [2m — n — l,n + l — m] (34) 

Substitution of this result into (^) immediately yields the result V/C„(y) = /C„+i(y) and so 
completes our proof. 
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